Hi again! We hope you are excited to see how the concepts you just learned about work in practice! Let’s dive right in and go through some examples!

Setup

install packages

If you’re running Mac or Linux, the previous command to install sf may not work first time. These operating systems (OSs) have ‘systems requirements’ that are described in the package’s README. Various OS-specific instructions can be found online, such as the article Installation of R 4.0 on Ubuntu 20.04 on the blog rtask.thinkr.fr.

#install.packages("sf") 
#install.packages("spData")

Load packages

sf - Simple Features for R

spData - Diverse spatial datasets for demonstrating, benchmarking and teaching spatial data analysis. It includes R data of class sf (defined by the package ‘sf’)

library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(ggplot2)
library(patchwork)

Wold data

The object loaded is a sf object containing a world map data from Natural Earth with a few variables from World Bank

# Alternative 1
world = spData::world

# Alternative 2
world = read_sf(system.file("shapes/world.gpkg", package = "spData"))
colnames(world)
##  [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"     
##  [7] "area_km2"  "pop"       "lifeExp"   "gdpPercap" "geom"
head(world$name_long, 20)
##  [1] "Fiji"                             "Tanzania"                        
##  [3] "Western Sahara"                   "Canada"                          
##  [5] "United States"                    "Kazakhstan"                      
##  [7] "Uzbekistan"                       "Papua New Guinea"                
##  [9] "Indonesia"                        "Argentina"                       
## [11] "Chile"                            "Democratic Republic of the Congo"
## [13] "Somalia"                          "Kenya"                           
## [15] "Sudan"                            "Chad"                            
## [17] "Haiti"                            "Dominican Republic"              
## [19] "Russian Federation"               "Bahamas"

CRS with SF

methods in SF

methods(class = "sf")
##  [1] [                     [[<-                  $<-                  
##  [4] aggregate             as.data.frame         cbind                
##  [7] coerce                dbDataType            dbWriteTable         
## [10] filter                identify              initialize           
## [13] merge                 plot                  print                
## [16] rbind                 show                  slotsFromS3          
## [19] st_agr                st_agr<-              st_area              
## [22] st_as_s2              st_as_sf              st_bbox              
## [25] st_boundary           st_buffer             st_cast              
## [28] st_centroid           st_collection_extract st_convex_hull       
## [31] st_coordinates        st_crop               st_crs               
## [34] st_crs<-              st_difference         st_filter            
## [37] st_geometry           st_geometry<-         st_inscribed_circle  
## [40] st_interpolate_aw     st_intersection       st_intersects        
## [43] st_is_valid           st_is                 st_join              
## [46] st_line_merge         st_m_range            st_make_valid        
## [49] st_nearest_points     st_node               st_normalize         
## [52] st_point_on_surface   st_polygonize         st_precision         
## [55] st_reverse            st_sample             st_segmentize        
## [58] st_set_precision      st_shift_longitude    st_simplify          
## [61] st_snap               st_sym_difference     st_transform         
## [64] st_triangulate        st_union              st_voronoi           
## [67] st_wrap_dateline      st_write              st_z_range           
## [70] st_zm                 transform            
## see '?methods' for accessing help and source code

Retrieve CRS

st_crs(world)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

Retrieve additional information about CRS

st_crs(world)$IsGeographic
## [1] TRUE
st_crs(world)$units_gdal 
## [1] "degree"
st_crs(world)$srid 
## [1] "EPSG:4326"
st_crs(world)$proj4string 
## [1] "+proj=longlat +datum=WGS84 +no_defs"

Set CRS

world_wgs84 = st_set_crs(world, "EPSG:4326") # set CRS to WGS84
st_crs(world_wgs84)$srid 
## [1] "EPSG:4326"

Projections

simple plotting with base R

Here we have a geographic coordinate system

#plot(world)
plot(world['continent'])

Projections with sf

head(sf_proj_info(type = "proj"))
##         name                  description
## 1 adams_hemi Adams Hemisphere in a Square
## 2  adams_ws1    Adams World in a Square I
## 3  adams_ws2   Adams World in a Square II
## 4        aea            Albers Equal Area
## 5       aeqd        Azimuthal Equidistant
## 6     affine        Affine transformation

Here we have two projected coordinate systems

# Robinson Projections
world_robin <- st_transform(world, "+proj=robin")
plot(world_robin['lifeExp']) 

# Mollweide Projections
world_moll <- st_transform(world, "+proj=moll")
plot(world_moll['gdpPercap'])

Here we can the the see the distortions in each of the systems for the USA

plot(world_robin['geom'][world_robin$name_long == 'United States', ])

plot(world_moll['geom'][world_moll$name_long == 'United States', ])

plot(world['geom'][world$name_long == 'United States', ])

Finding the right coordinate systems

We recommend to check out https://epsg.io for details about coordinate systems

Here we compare the US National Atlas Equal Area projection with WGS 84 coordinate systems

# the package data is stored in NAD83
st_crs(us_states)$srid 
## [1] "EPSG:4269"
st_crs(us_states)
## Coordinate Reference System:
##   User input: EPSG:4269 
##   wkt:
## GEOGCS["NAD83",
##     DATUM["North_American_Datum_1983",
##         SPHEROID["GRS 1980",6378137,298.257222101,
##             AUTHORITY["EPSG","7019"]],
##         TOWGS84[0,0,0,0,0,0,0],
##         AUTHORITY["EPSG","6269"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4269"]]
# EPSG:2163 - US National Atlas Equal Area
p1 <- ggplot(data = us_states) +
    geom_sf() +
    coord_sf(crs = st_crs(2163)) + theme_bw()

# EPSG:4326 - WGS 84 
p2 <- ggplot(data = us_states) +
    geom_sf() +
    coord_sf(crs = st_crs(4326)) + theme_bw()

library(patchwork)
p1 + p2

Test EPSG:2163 on world map

ggplot(data = world) +
    geom_sf() +
    coord_sf(crs = st_crs(2163)) + theme_bw()

Working with spatical data and plots

Australia vs. USA

Here we are combining two plots of Australia and the USA

us_sf <- world['geom'][world$name_long == 'United States', ]
us_sf <- st_transform(us_sf, "EPSG:4326")

aus_sf <- world['geom'][world$name_long == 'Australia', ]
aus_sf <- st_transform(aus_sf, "EPSG:4326")

us <- ggplot() +
  geom_sf(data = us_sf)

aus <- ggplot() +
  geom_sf(data = aus_sf)


us + aus

st_crs(us_sf)$srid 
## [1] "EPSG:4326"
st_crs(aus_sf)$srid 
## [1] "EPSG:4326"

Now we are plotting Australia and the USA in a single plot.

us_sf <- st_transform(us_sf, 'EPSG:4326')
aus_sf <- st_transform(aus_sf, 'EPSG:4326')

st_crs(us_sf)$units_gdal 
## [1] "degree"
ggplot()+
  geom_sf(data = us_sf)+
  geom_sf(data = aus_sf)+ 
  coord_sf(xlim = c(-130, 150), ylim = c(-50, 50), expand = FALSE) +
  theme_void()

We can move Australia next to the USA. Keep in mind that we are still using WGS80 / EPSG:4326. That means we have distortions in our plot.

# copy the object
map_au_moved <- aus_sf
# transform (move) the object by -180 degrees + 60 degrees (north/west)
st_geometry(map_au_moved) <- st_geometry(map_au_moved) + c(-180, 60)
# copy the coordinate reference systems
st_crs(map_au_moved) <- st_crs(aus_sf)

ggplot()+
  geom_sf(data = us_sf)+
  geom_sf(data = map_au_moved)+ 
  coord_sf(xlim = c(-130, -20), ylim = c(15, 50), expand = FALSE) +
  theme_void()
## Warning in st_cast.GEOMETRYCOLLECTION(X[[i]], ...): only first part of
## geometrycollection is retained

To minimize distortion we can use local datums. Then we use these projections in coordinate system that reflect the same scale

us_sf <- st_transform(us_sf, 'EPSG:2163')   # EPSG:2163 - US National Atlas Equal Area
aus_sf <- st_transform(aus_sf, 'EPSG:3112') # EPSG:3112 - GDA94 / Geoscience Australia Lambert

st_crs(us_sf)$units_gdal 
## [1] "metre"
st_crs(aus_sf)$units_gdal 
## [1] "metre"
ggplot()+
  geom_sf(data = us_sf)+
  # coordinate system with 6.000 km width and 4.000 km height
  coord_sf(xlim = c(-3000000,3000000), ylim = c(-3000000,1000000)) +
  theme_bw()

ggplot()+
  geom_sf(data = aus_sf)+
  # coordinate system with 6.000 km width and 4.000 km height
  coord_sf(xlim = c(-3000000,3000000), ylim = c(-5000000,-1000000)) +
  theme_bw()

Displaying and cropping spatial dataset

Without modifying the spatial dataset

Displaying a map of Europe by specifying the longitudes and latitudes with coord_sf() Notice: coord_sf() is applied on the whole world data.

ggplot() + geom_sf(data = world) +
    coord_sf(xlim = c(-25, 50), ylim = c(30, 73), expand = FALSE) + theme_bw()

With modifying the spatial dataset

Here we filter our World dataset with ‘Europe’ The European overseas territories and Russia are not exactly what we wanted.

europe <-world['geom'][world$continent == 'Europe', ]
ggplot() + geom_sf(data = europe) + theme_bw()

Now we use the cropping function st_crop() to cut our dataset.

europe_cropped <- st_crop(europe, xmin = -25, xmax = 50,
                                  ymin =  30, ymax = 73)

ggplot() + geom_sf(data = europe_cropped) + theme_bw()

Locate Hertie School and create our own geometries

To locate the Hertie School on our map we first have to retrieve the coordinates. Then we can create a point geometry and plot it on top of our map.

coord_hertie <- c(13.389225005174465, 52.512803179279636) 

hertie_school = st_sfc(st_point(coord_hertie), crs = 4326)

ggplot() + geom_sf(data = europe_cropped) +
           geom_sf(data = hertie_school, color = 'red') + 
           theme_bw()

With st_buffer we can increase our geometry

hertie_buff = st_buffer(hertie_school, dist = 50 * 1000)

ggplot() + geom_sf(data = europe_cropped) +
           geom_sf(data = hertie_buff, fill = 'red', alpha = 0.5, color = 'red') + 
           theme_bw()

The st_buffer also allows us to draw wider circles around our buffer. Since our CRS is still WGS80 we have distortions in Europe. Therefore, the circles are not perfect.

buff_500km  = st_buffer(hertie_school, dist = 500 * 1000)
buff_1000km = st_buffer(hertie_school, dist = 1000 * 1000)
buff_1500km = st_buffer(hertie_school, dist = 1500 * 1000)

ggplot() + geom_sf(data = europe_cropped) +
           geom_sf(data = hertie_buff, fill = 'red', alpha = 0.5, color = 'red') +   
           geom_sf(data = buff_500km, fill = 'grey', alpha = 0.1) +
           geom_sf(data = buff_1000km, fill = 'grey', alpha = 0.1, ) + 
           geom_sf(data = buff_1500km, fill = 'grey', alpha = 0.1, ) + 
           theme_bw()

st_crs(europe_cropped)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

To archive better circles we have to use a datum that is more suitable for Europe. For examples EPSG:5643 is recommended for central Europe.

new_europe <- st_transform(europe_cropped, 'EPSG:5643') # set CRS
new_buff <- st_transform(hertie_buff, 'EPSG:5643')      # set CRS
buff_500km  <- st_transform(buff_500km, 'EPSG:5643')    # set CRS
buff_1000km <- st_transform(buff_1000km, 'EPSG:5643')   # set CRS
buff_1500km <- st_transform(buff_1500km, 'EPSG:5643')   # set CRS




ggplot() + geom_sf(data = new_europe) +
           geom_sf(data = new_buff, fill = 'red', alpha = 0.5, color = 'red') +   
           geom_sf(data = buff_500km, fill = 'grey', alpha = 0.1) +
           geom_sf(data = buff_1000km, fill = 'grey', alpha = 0.1, ) + 
           geom_sf(data = buff_1500km, fill = 'grey', alpha = 0.1, ) + 
           theme_bw()

st_crs(new_europe)$units_gdal 
## [1] "metre"
st_crs(new_buff)$units_gdal 
## [1] "metre"
st_crs(new_europe)$IsGeographic
## [1] FALSE
st_crs(new_buff)$IsGeographic
## [1] FALSE

Excercise 1

Try to retrive the spatial data for Canada from World country polygons (from the spData package)

Excercise 2

Try to find a suitable projection or datum for Canada on https://epsg.io and use it for plotting Canada in ggplot2.

Excercise 3

Try to find the coordinates for Canadas capital (Ottawa) and create a point geometry. Plot the point to that map you presivouly created.

Hopefully you learned something new! Feel free to contact us in case you have further questions and enjoy the rest of our interesting workshop day!

Sources